
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%       Test: Choice of Specification for the Bidders' Private Values     %
%                      Linear versus Exponential                          %
%                        Monte Carlo Experiment                           %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all

L2=200; % number of auctions with N=2
L3=200; % number of auctions with N=3
L=L2+L3;
N=[2,3];
alpha=0.50:0.10:0.80; % quantile levels alpha
psi2=ones(L2,1)*(N(1)*(alpha.^(N(1)-1))-(N(1)-1)*(alpha.^N(1))); % Psi function given N=2
psi3=ones(L3,1)*(N(2)*(alpha.^(N(2)-1))-(N(2)-1)*(alpha.^N(2))); % Psi function given N=3
psi = [psi2;psi3];
m=size(alpha,2);

% Auxiliary Parameters
nsim=10000; % Number of simulations
B=500; % Bootstrap replications
toler=0.1; % Parameter of tolerance for estimation
lambda_star=0.5846; % Parameter that satifies the equivalence property between a linear and an exponential specification (obtained in a previous simulation exercise - Code: Vuong_lambda_star)
Rej1=0; % Auxiliary sum for rejection probability
Rej5=0;
Rej10=0;

gammai=[0,0.5]'; % Initial guess for coefficients 

% Simulation 
for i=1:nsim
    
   % Data Generating Process
   eps2=unifrnd(0,1,L2*N(1),1); % Regression Errors N=2
   eps2=reshape(eps2,N(1),L2);
   eps2=sort(eps2,'descend');
   eps2=eps2(2,:);
   eps3=unifrnd(0,1,L3*N(2),1); % Regression Errors N=3
   eps3=reshape(eps3,N(2),L3);
   eps3=sort(eps3,'descend');
   eps3=eps3(2,:);
   eps=[eps2';eps3'];
   
  
   x2=[ones(L2,1),normrnd(0,1,L2,1)]; % Covariates N=2
   x3=[ones(L3,1),normrnd(0,1,L3,1)]; % Covariates N=3
   x=[x2;x3];
   p=size(x,2);
   
   wbL=x*gammai+eps; % Winning bid specification under a linear model for private values
   wbE=exp(x*gammai+eps); % Winning bid specification under a exponential model for private values
   wb=lambda_star*wbL+(1-lambda_star)*wbE; % Linear combination of both specifications
   dat=[wb,x]; % Matrix with data
   
   % Computation of the choice of specification test statistic 
    for k=1:m
           % objective function under H0
           [objE ~]=MM_Exp(psi(:,k), p, wb, x, toler, gammai);
           % Objective function under H1
           [objL ~]=MM(psi(:,k), p, wb, x, toler, gammai);
           Vuong(k,:)=objE-objL; 
    end

   Vuongmax=sqrt(L)*max(Vuong); % test statistic
     
% Pairwise Bootstrap of the choice of specification test statistic
for b=1:B
    datb2 = dat(randsample(1:L2,L2,true),:); % Resampling of the matrix of data in the subsample for N=2
    datb3 = dat(randsample(L2+1:L,L3,true),:); % Resampling of the matrix of data in the subsample for N=3
    wb2 = datb2(:,1); % Resampled winning bids N=2
    wb3 = datb3(:,1); % Resampled winning bids N=2
    wb=[wb2;wb3];
    x2 = datb2(:,[2 3]); % Resampled covariate N=3
    x3 = datb3(:,[2 3]); % Resampled covariate N=3   
    x=[x2;x3];
    
    for k=1:m
        % objective function under H0
        [objEb ~]=MM_Exp(psi(:,k), p, wb, x, toler, gammai);
        % objective function under H1
        [objLb ~]=MM(psi(:,k), p, wb, x, toler, gammai);
        Bootvuong(k,:)=sqrt(L)*((objEb-objLb)-Vuong(k,:));
    end
   
   Bootvuongmax(:,b)=max(Bootvuong); % bootstraped test statistic
            
end
    % Critical Values
    c99=quantile(Bootvuongmax',0.99);
    c95=quantile(Bootvuongmax',0.95);
    c90=quantile(Bootvuongmax',0.90);

    % Computation of rejection probability 
    if (Vuongmax>c99)
        Rej1=Rej1+1;
    end
    if (Vuongmax>c95)
        Rej5=Rej5+1;
    end
    rej5=Rej5/nsim;

    if (Vuongmax>c90)
        Rej10=Rej10+1;
    end

    i
end

% Rejection probability
Rej1=Rej1/nsim
Rej5=Rej5/nsim
Rej10=Rej10/nsim